home *** CD-ROM | disk | FTP | other *** search
/ Disc to the Future 2 / Disc to the Future Part II Programmer's Reference (Wayzata Technology)(6013)(1992).bin / MAC / THINKC / 4_0 / VIVIDUS / QD3D.SIT / qd3d / geometry.c < prev    next >
Text File  |  1991-10-09  |  5KB  |  252 lines
  1. /*    ======================================================================
  2.  
  3.     Supplemental geometry routines.
  4.     
  5.     This is part of the qd3d Vividus Source Code Library.  See the
  6.     extern qd3d.doc documentation file for usage.  See individual
  7.     routines for routine documentation.
  8.     
  9.     Copyright 1991 by Vividus Consulting.
  10.     
  11.     This is not public domain source code.  You may not copy and
  12.     paste from this source code.  Read your Vividus Licensing
  13.     agreement for details and other restrictions.
  14.     
  15.     Note:    At present, these routines are only supported as needed by
  16.             hedra.
  17.  
  18.     ======================================================================    */
  19.  
  20. #include    <math.h>
  21. #include    <vect.h>
  22. #include    <vect.g>
  23. #include    "geometry.h"
  24.  
  25. static double    twoPI = 6.283185307179586477;
  26. static double    halfPI = 1.570796326794896619;
  27.  
  28. static double    two = 2.0;
  29. static double    four = 4.0;
  30. static double    zerop = .000001;    /*    Zero on the positive side    */
  31.  
  32. static double
  33. DizzyAngle(
  34.     vector    *a,
  35.     int        plane);
  36.  
  37. static double
  38. DizzyDelta(
  39.     double    b,
  40.     double    a);
  41.  
  42. /*    ======================================================================    */
  43.  
  44. #define    yzPlane    0
  45. #define    xzPlane    1
  46. #define    xyPlane 2
  47.  
  48. Boolean
  49. PtInPoly(
  50.     vector    *pt,        /*    Point on poly plane to see if in poly.        */
  51.     int        n,            /*    Number of verts in poly.                    */
  52.     vector    verts[],    /*    Vertices of polygon.                        */
  53.     vector    *norm)        /*    Normal of polygon.                            */
  54. /*
  55.     Returns true if the point (presumably already on the plane of the
  56.     polygon is in the polygon.
  57.     
  58.     Note:
  59.         This only works if the polygon is planar.
  60.         
  61.         The intersection of the ray with the poly plane must already
  62.         be found.
  63. */
  64. {
  65.     vector    a;
  66.     double    theta = 0.0;
  67.     double    thetasum = 0.0;
  68.     double    thetanew;
  69.     double    theta0;
  70.     int        i;
  71.     int        orientation;
  72.     
  73.     if    (norm->z > norm->y) {
  74.         if (norm->z > norm->x)
  75.             orientation = xyPlane;
  76.         else
  77.             orientation = yzPlane;
  78.     } else {
  79.         if (norm->y > norm->x)
  80.             orientation = xzPlane;
  81.         else
  82.             orientation = yzPlane;
  83.     }
  84.     
  85.     vvect(pt, &verts[0], &a);
  86.     theta = theta0 = DizzyAngle(&a, orientation);
  87.     
  88.     for (i = 1; i < n; i++) {
  89.         vvect( pt, &verts[i], &a);
  90.         thetanew = DizzyAngle(&a, orientation);
  91.         thetasum += DizzyDelta(thetanew, theta);
  92.         theta = thetanew;
  93.     }
  94.     thetasum += DizzyDelta(theta0, theta);
  95.             
  96.     if ( (thetasum > PI) || (thetasum < -PI) )
  97.         return(true);
  98.     else
  99.         return(false);
  100. }
  101.  
  102. Boolean
  103. RayXPlane(
  104.     Ray        *r,
  105.     Plane    *p,
  106.     double    *t,
  107.     vector    *pt)
  108. /*
  109.     Returns true if the ray intersects the plane and the point and paramter of
  110.     the intersections.
  111. */
  112. {
  113.     vector    v;
  114.     
  115.     vsub(&p->pt, &r->b, &v);
  116.     *t = vdot(&r->m, &p->normal);
  117.     if (*t == 0.0)
  118.         return(false);
  119.         
  120.     *t = vdot(&v, &p->normal)/(*t);
  121.     
  122.     vscale(*t, &r->m, &v);
  123.     vadd(&v, &r->b, pt);
  124.     
  125.     return(true);
  126. }
  127.  
  128. int
  129. RayXEllipsoid(
  130.     Ray            *r,
  131.     Ellipsoid    *e,
  132.     double        *t)
  133. /*
  134.     Returns false if there isn't an intersection.  If there is, the number of
  135.     intersections is returned; ed, eb, &t reflect the point, normal, & t of the ray.
  136.     
  137.     Note:        Intersections can only occur with positive values of t.
  138.     Warning:    Presently the only ellipsoid supported is a sphere.
  139. */
  140. {
  141. /*    double    a = r->m.x * r->m.x + r->m.y * r->m.y + r->m.z * r->m.z;    */
  142.     double    a = 1.0;    /*    rm is a unit vector.    */
  143.     double    b = two * r->m.x * (r->b.x - e->center.x) +
  144.                 two * r->m.y * (r->b.y - e->center.y) +
  145.                 two * r->m.z * (r->b.z - e->center.z);
  146.     double    c = e->center.x * e->center.x +
  147.                 e->center.y * e->center.y +
  148.                 e->center.z * e->center.z +
  149.                 r->b.x * r->b.x + r->b.y * r->b.y + r->b.z * r->b.z -
  150.                 two * (    e->center.x * r->b.x +
  151.                         e->center.y * r->b.y +
  152.                         e->center.z * r->b.z) -
  153.                 e->dimensions.x * e->dimensions.x;
  154.     double    d = b * b - four * a * c;
  155.     double    t1, t2;
  156.     vector    p1, p2;
  157.     int        retval = 0;
  158.     
  159.     
  160.     if (d < 0.0)
  161.         return (retval);
  162.     
  163.     if (d == 0.0) {
  164.         t1 = -b / (two * a);
  165.         if (t1 > zerop)
  166.             retval = 1;
  167.     } else {
  168.         t1 = (- b - sqrt(d))/(two * a);
  169.         t2 = (- b + sqrt(d))/(two * a);
  170.         if (t1 > zerop)
  171.             retval = 2;
  172.         else {
  173.             if (t2 > zerop)
  174.                 retval = 1;
  175.             else
  176.                 return(0);
  177.         }
  178.     }
  179.     
  180.     /*    t1 (if it is +) must be the closer root    */
  181.     if (retval == 1)
  182.         *t = t2;        /*    t1 was negative.    */
  183.     else
  184.         *t = t1;
  185.     return (retval);
  186. }
  187.  
  188. void
  189. EllipsoidFromBox(Bound3dBox *bb, Ellipsoid *e)
  190. /*
  191.     Finds the smallest ellipsoid containing the box.  The ellipsoid is
  192.     centered about the box.
  193.     
  194.     Note:    Presently only spherical ellipsoids are created.
  195. */
  196. {
  197.     vector    a;
  198.     double    s = 1.0/2.0;
  199.     
  200.     vadd(&bb->max, &bb->min, &a);
  201.     vscale(s, &a, &e->center);
  202.     s = vdist(&e->center, &bb->max);
  203.     vscale(s, &vone, &e->dimensions);
  204. }
  205.  
  206. /*    ======================================================================    */
  207.  
  208. static double
  209. DizzyAngle(
  210.     vector    *a,
  211.     int        plane)
  212. /*
  213.     Returns the dizzy angle.
  214. */
  215. {
  216.     switch(plane) {
  217.         case yzPlane:
  218.             return(atan2(a->y, a->z));
  219.         case xzPlane:
  220.             return(atan2(a->x, a->z));
  221.         case xyPlane:
  222.             return(atan2(a->x, a->y));
  223.     }
  224. }
  225.  
  226. static double
  227. DizzyDelta(
  228.     double    b,
  229.     double    a)
  230. /*
  231.     deltat theta from angle a to angle b.
  232.     
  233.     Note this returns a value between -╣ and ╣.  It is necessary for when the vertices
  234.     "wrap around" the -╣, ╣ boundary.
  235. */
  236. {
  237.     double    dt = b - a;
  238.     
  239. #if    1    /*    Should be turned off once we're sure -╣ <= atan2 <= ╣    */
  240.     if ((a > PI) || (a < -PI) || (b < -PI) || (b > PI))
  241.         DebugStr("\pdizzydelta domain error");
  242. #endif
  243.  
  244.     if (dt > PI)
  245.         dt -= twoPI;
  246.     if (dt < -PI)
  247.         dt += twoPI;
  248.     
  249.     return(dt);
  250. }
  251.  
  252.